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Abstract 

Solving the 4-d Einstein equations as evolution in time requires solving equations of two types: 
the four elliptic initial data (constraint) equations, followed by the six second order evolution equa- 
tions. Analytically the constraint equations remain solved under the action of the evolution, and 
one approach is to simply monitor them (unconstrained evolution). Since computational solution of 
differential equations introduces almost inevitable errors, it is clearly "more correct" to introduce 
a scheme which actively maintains the constraints by solution (constrained evolution). This has 
shown promise in computational settings, but the analysis of the resulting mixed elliptic hyper- 
bolic method has not been completely carried out. We present such an analysis for one method of 
constrained evolution, applied to a simple vacuum system, linearized gravitational waves. 

We begin with a study of the hyperbolicity of the unconstrained Einstein equations. (Because the 
study of hyperbolicity deals only with the highest derivative order in the equations, linearization 
loses no essential details.) We then give explicit analytical construction of the effect of initial 
data setting and constrained evolution for linearized gravitational waves. While this is clearly a 
toy model with regard to constrained evolution, certain interesting features are found which have 
relevance to the full nonlinear Einstein equations. 

PACS numbers: 
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I. INTRODUCTION 



Binary black hole systems are expected to be the strongest possible astrophysical gravi- 
tational wave sources. In the final moments of stellar mass black hole inspiral, the radiation 
will be detectable in the current (LIGO-class) detectors. If the total binary mass is of the 
order of 1OM , the moment of final plunge to coalescence will emit a signal detectable by the 
current generation of detectors from very distant (Gpc) sources. The merger of supermas- 
sive black holes in the center of galaxies will be the very dominant signal in the spaceborne 
LISA detector, and detectable out to large redshift. Simulation of these mergers will play an 
important part in the prediction, detection, and the analysis of their gravitational signals in 
gravitational wave detectors. To do so requires a correct formalism which does not generate 
spurious singularities during the attempted simulation. 

Recent work at Texas [J P] has found that constrained 3-d evolution can produce sub- 
stantially stabilized long-term single black hole simulations, stimulating interest in con- 
strained evolution. However, it is true that to date there has been no analysis that addresses 
the behavior of constrained evolution for computational relativity. That is the purpose of 
this paper. (Note that unlike the situation in electromagnetism, where a discrete simplectic 
calculus can be written which assures conservation of a discretized version of the E&M con- 
straints, there appears to be no equivalent formulation for general relativity that conserves 
discretized versions of the gravitational constraints jjj.) 



II. 3 + 1 FORMULATION OF EINSTEIN EQUATIONS 

We take a Cauchy formulation (3+1) of the ADM type, after Arnowitt, Deser, and 
Misner Q| . In such a method the 3- metric and its momentum (the extrinsic curvature) 
are specified at one initial time on a spacelike hypersurface, and evolved into the future. The 
ADM metric is 

ds 2 = -(a 2 - fop) dt 2 + 2pi dt dx l + 9ij dx i dx j (1) 

where a is the lapse function and /3 l is the shift 3- vector; functions that relate the coordinates 
on each hypersurface to each other. [16] 

[16] Latin indices run 1, 2, 3 and are lowered and raised by <jry and its 3-d inverse g lJ . 
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The Einstein field equations contain both hyperbolic evolution equations, and elliptic 
constraint equations. The constraint equations for vacuum in the ADM decomposition are: 

H=±[R-K ij K i > + K 2 ] = 0, (2) 

H* = Vj (K ij - g lj K) = 0. (3) 

Here R is the 3-d Ricci scalar constructed from the 3-metric, and Vj is the 3-d covariant 
derivative compatible with g^. Initial data must satisfy these constraint equations; one may 
not freely specify all components of gij and K^. 

The evolution equations from the Einstein system are 

dtgtj = -2aKij + V^A + V^-, (4) 

and 

dtKij = -ViVja + a [Rtj - 2K ia K a j + KK i3 \ + (3 k K ijjk + K kj (3* + K ik p$ , (5) 

where Rij is the 3-d Ricci tensor. 

We refer to this form of the Einstein equations as of ADM type, referring to the funda- 
mental development this specific form is called the g - K form. Here, Eq. (j2J)-Eq. Q, the 
constraint equations, are the vacuum Einstein equations 4 G o = and 4 G i = respectively, 
for a 4- metric of form Eq. (0). Eq. (j3J)-Eq. (jSJ), the evolution equations, are a first order 
form of the vacuum Einstein equations A R%j = 0. The methods decribed here can be easily 
extended to linearized versions of some other formulations of Einstein's equations, and one 
of these is explicitly addressed below. 

III. DATA FORM 

In this paper we consider only linearized gravitational waves (but see Appendix A). 

We write the spatial metric as g a b = S a b + h ab . We work with zero shift vector: (3 l = 
(cf Q]), but we consider a lapse a which may differ from the flat background value by a 
function of first order in the linearization: a = 1 + a\+0{h 2 ab ) . In particular we choose 
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below a densitized lapse: 



a = g b 



= l + b(h xx + h yy + h zz ) + 0(h 2 ab ), (6) 

where the quantity g is the determinant of the 3— metric, and b is a chosen constant. Gen- 
erally, densitized lapse can involve multiplication by a fixed function of the coordinates, but 
we assume linearization from flat space, so this function here must be unity. 
For our linearized case, Eq. (HJ) reads 

d t h l3 = -2K ijy (7) 

and Eq. (j^J) is: 

d t Kij = = -^^(didmhij + didjh lm - didih mj - djdih m ^ - didja. (8) 

Eq. © can, with Eq. (|7jl . be read as a second-order in time equation for hij\ 

dthij = a 2 5 lm (did m h i:j + didjh lm - didih mj - djdih m ^) + 2<9 i <9 j a, (9) 

We further assume a 1— dimensional system, so that all quantities are functions of 
(x 1 ,^) = (x, t), where a; is a spatial variable. It is then profitable to write out the indi- 
vidual components explicitly. We shall find that the different components xx (longitudinal- 
longitudinal: LL); xy & xz (transverse- longitudinal: TL); yy + zz (transverse trace); and 
yy — zz & yz (transverse-traceless: TT) have, as expected, different behavior. This explicit 
formulation parallels, in a less sophisticated but more transparent manner, the decomposi- 
tion in 1^ | . 

The second-order evolution equations, written explicitly in terms of these variables are: 

<%h xx = d x 2 (h yy + h zz ) + 2didja (10) 

d 2 h xy = (11) 

d 2 t h xz = (12) 

d 2 (h yy + h zz ) = d x 2 (h yy + h zz ) (13) 
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d 2 t (hyy-h ZZ ) = d X 2 (hyy-h ZZ ) (U) 

d 2 t {h yz ) = d x \hy z ). (15) 

The constraint equations become: 

H = -d 2 x (h yy + h zz ), (16) 

H x = -2d x (Kyy + K zz ), (17) 

Hy = 2d X K X y, (18) 

H z = 2d x K xz . (19) 

Note that since the spatial metric is 5 a b, we drop the distinction between raised and 
lowered indices in these explicit Equations (jlUj) - (|19j) . We relate the constraint equations 
ffT7f) - (HHD to d t h i: j via Eq. J7J) above. 

IV. INITIAL DATA SETTING AND CONSTRAINED EVOLUTION 

We first give the full, nonlinear data setting (constraint solving) procedure, then specialize 
in the next section to our linearized plane wave case. 

We adopt the confonnal t_e-t r ace,e. m etKod of Yo rk and collaborator Q-Q 
which consists of a conformal decomposition with a scalar <fi that adjusts the Hamiltonian 
constraint, and a vector potential w l that adjusts the longitudinal components of the extrin- 
sic curvature. The constraint equations are then solved for these new quantities such that 
the complete solution fully satisfies the constraints. 

For initial data setting we pose a trial metric and a trial trace-subtracted extrinsic cur- 
vature taken as conformal trial functions and A l K 

The physical metric, g^-, and the trace- free part of the extrinsic curvature, Ay, are related 
to the background fields through a conformal factor 

9ij = (j> 4 9ij, (20) 
A ij = (/r 10 (i lJ + {Iwf), (21) 



where is the conformal factor, and (Iw) will be used to cancel any possible longitudinal 
contribution in A. w % is a vector potential, and 

(Iwf = VV + V j w i - -g ij V k w k . (22) 

(Here is the covariant derivative in the conformal background space.) The trace K is 
not corrected: 

K = K. (23) 

Writing the full, nonlinear, Hamiltonian and momentum constraint equations in terms of 
the quantities in Eqs. (|2~T|) - (J23j) . we obtain four coupled elliptic equations for the fields (p 
and w l P|: 

V 2 = (1/8) W +\k 2 4? ~ 

+ (fwfXAij + {lw)ij)), (24) 

S/^lwf = -tf'ttPVjK -VjA ij . (25) 

After the data have been set, the evolution begins. A well designed evolution scheme 
evolves to the next timestep. In an unconstrained evolution, this completes the time step. 
In our constrained scheme, the metric and extrinsic curvature so determined are taken as 
intermediate values in the middle of a timestep and thus as conformal trial functions gij and 
A*. 

The constraint Equations (|24j ) .(|25 p are solved with these trial functions to complete each 
time update step. The resulting solved gij and Kij are taken as the data for the next time 
update. Notice that these equations require no specific gauge choice; the constraint solution 
is independent of the gauge functions a and (3 l . A similar approach also can be applied to 
other formulations which generally have a larger number of constraints. 

V. INITIAL DATA SETTING FOR THE LINEARIZED SYSTEM 

We assume a choice of h a b and K a b. We then solve the linearized version of the constraint 
Equations (|24j) . f!25[) . obtaining and w l . As in Section III we assume these variables are 
functions only of x (not of time, because this is the intial data problem, which is solved at 
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one instant of time). Then since the extrinsic curvature vanishes in the (flat) background, 
the linearized version of the Hamiltonian constraint (J2*lj) is simply 



%4> = \r, 

= -\dl{h yy + h zz ). (26) 

The linearized version of the momentum constraints (|25|). written explicitly for each com- 
ponent, are: 

y| = \K X - A™ x . (27) 

w*>* x = -Ay%. (28) 

w z '% = -A*%. (29) 
A. Elliptic Equation Boundary Conditions 

A solution of the elliptic constraint equations requires that boundary data be specified. 
For this demonstration of the technique, analytic integration in our 1-dimensional system 
will lead to arbitrary integration constants. Specific choices for these constants is equivalent 
to setting 6 and w t,x to specific boundary values; this will be amplified in the Appendix, 
where we give a brief discussion of the effect of such boundary condition choices. For 
instance, we will see that setting the integration constants to zero leads to simplest analytical 
results. However this is somewhat at variance with the process we use in full 3-dimensional 
simulations to handle black holes, where we choose conditions 6 = 1 and w l = at the inner 



boundaries (the excision boundaries of the black holes); and mixed (Robin) ll| conditions 
at the outer boundaries for black holes 

B. Constraint Solution for Linear Waves 

In the linearized regime the Equations (}2TH) -(f2"9" j) have decoupled. Each of these equations 
is an ordinary differential equation, and is straightforwardly integrated. We obtain solutions: 



'i = -\{h V y + hzz) + Ax + B. (30) 
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\w x > x = \k-A x * + C X) (31) 
w y,x = _Av* + C y , (32) 
w z,x = _fr* + c ^ (33) 

where A,B,Ci are integration constants. Note that since w l is a potential, we need only its 
first derivatives. 

^From Eq. (|2Tj) . the effect of the Hamiltonian constraint solution for <fi is: 

9ab = 3ab + h ab 

= exp40[5 ab + ha/,] 

= (1 + 40 1 )[5 o6 + M + 0(/i ab 2 ). (34) 

At the linearized level this conformal factor does not affect the off-diagonal components, 
so we obtain: 

Kx = h xx - ^(h yy + h zz ) + 4Ax + 4B, (35) 

h X z = h xz , (37) 
h yy + h zz = + 8Ax + 8B, (38) 

foyy hzz h"yy h Z zi (39) 

h yz = h yz - (40) 

Notice that the TT (yy — zz & yz) components are unchanged, but the transverse trace 
(yy + zz) is set to zero modulo integration constants, with a similar subtraction of the LL 
(xx) component, which, however, generally leaves h xx nonzero. With only this one variable 
it is not surprising that the TL terms (zx & yx) are left unchanged. 

The effect of solution of the momentum constraints is to modify the extrinsic curvature 
by adding (£w) . The components w y and w z set A xy and A xz to zero, respectively (modulo 
integration constants). w x modifies all the diagonal components of Aif A xx = A xx + |u> x,x , 
A yy = Ayy — |iy a ' ,x , A zz = A zz — ^w x ,x ; note that A ab is traceless since A ab is. The results: 
A xx = \K, Ayy = Ay y + \A XX - \ K , A zz = A zz + \A XX - \ K. Written in terms of K ab rather 
than A ab , and restoring the integration constants, we find: 
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K 


= K + C x , 


(41) 


K 


— U Uy, 




K 


= o + c z , 


(43) 


yy 


= o - c x , 


(44) 


K — K 

yy Iy zz 


— k — k 

yy lv zzi 


(45) 


Kyz 




(46) 



Note, as assumed, that the trace, K, of the extrinsic curvature is not modified in the 
conformal constraint solution. 

By inserting these results fj33 J) - (f4T)J) into Equations (fTKjl - ffiTjl) . we can verify that the 
constraints are satisfied. 



VI. FREE EVOLUTION AND HYPERBOLICITY 

We now turn to the hyperbolicity of this system, assuming the data have been correctly 
set (the constraints are initially satisfied). Note that neither the LL metric component h xx 
nor its initial time derivative is generically zero. 



Kreiss and Ortiz |12| have shown that a second-order form, and its equivalent first or- 
der form have the same hyperbolicity classification. In pure second order form (no first 
derivatives) the classification is especially simple: If all the modes satisfy a wave equation 
with real nonzero propagation velocity, the system is strongly hyperbolic and well posed for 
Cauchy data. If any of the modes has zero propagation velocity, the system is at best weakly 
hyperbolic and is potentially ill posed. If any of the propagation velocities is imaginary the 
system is not hyperbolic, and is completely ill posed for Cauchy data. 

If a unit lapse is chosen so that d 2 x a = then the LL equation (fTU|) contains no spatial 
derivatives of h xx ; the propagation velocity for h xx is zero. Eq. fTT)|) then reflects at best 
weak hyperbolicity of the system; clearly it admits analytic linear growth in h xx , which can 
plainly lead to late-time difficulties. If, however, the lapse is densitized, with constant b > 
(cf Eq. ©), then Eq. (fTUj) becomes wavelike for h xx . If this were the only equation on the 
system this choice of lapse would have made the system strongly hyperbolic. 



However, one must investigate the entire set of evolution equations. Note that the TL 
components (fTTj) . (fT2|) of the field equations have no spatial derivatives, again signaling weak 
hyperbolicity, and analytic linear-growth solutions. Though the initial data can analytically 
set the slope to zero (if we choose C y = C z = 0), errors can arise in numerical evolution that 
triggers such linear growth. 

Thus the g — K version of these equations is weakly hyperbolic, and questions can arise 
as to the well-posedness of the solutions, even with densitized lapse. In refll3|, it is shown 

nn 

that a relatively simple extension that captures one of the main ideas of the BSSN|14I] |15| 
approach can cure the nonhyperbolicity of the TL components. This involves considering 
the combination 

fj = 5 kl d k h tj - U k h. (47) 

In that case one must treat the first order form of the equations because those for fj are 
first order in time. The evolution equation for fi is determined by differentiation of its 
definition. However the result is still a weakly hyperbolic system. By subtracting a multiple 
of the momentum constraint Hi to the right hand side of the equation for d t fi, a strongly 
hyperbolic system can finally be achieved. 

We will take a somewhat different approach here, as we discuss in the next section. The 
result just noted, that momentum constraint improves the hyperbolicity of the system, is of 
interest, because we will demonstrate momentum constraint solution, below. 



VII. CONSTRAINED EVOLUTION 



Constrained evolution as we described in Section IV above modifies the evolution. After 
each explicit step, the metric and extrinsic curvature (via Eq. (j7|l) are treated as conformal 
background quantities and are reset according to Eqs (jHo]) -(|4^ j) . Suppose in this Section that 
integration constants A, B, d are all zero. 

The transverse trace, which might have grown an error from the discretization of the 
equation, is reset to zero; cf Eq. (jHBJ)- Also, the momentum of this component is reset to 
zero (Eq. (jSSjl ). 

The h yy + h zz terms evolve to set the term h yy + h zz to zero in Eq. (J35|) . In this case the 
only terms on the right side of the d%h xx equation are the second derivatives of the lapse, 
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a. In setting data the trace of the extrinsic curvature, K^mai, is one of the data functions 
to be set, and with constrained data (and zero integration constants) K ini u a i gives d t h xx via 
Eq. 0. 

With densitized lapse, the h xx equation (J35j) is effectively a wave equation for h xx . Thus 
at least in this analytic discussion of the evolution, the linearized system appears completely 
controlled in a constrained, densitized-lapse evolution. 

If the lapse is chosen constant, but constraint solution Eq. ()38|) is imposed, then Eq. (fTUj) 
with Eq. ()41|) shows the the evolution is constant K (K = K xx is constant in time, though 
K x is not necessarily zero). While this causes no trouble in the linearized system, it will 
eventually violate the \h xx \ << 1 assumption. In the full nonlinear case this growth will 
eventually definitely cause problems. As we have noted, the nondensitized case is definitely 
weakly-hyp erb olic . 

Finally, notice that in this case that constraint solution sets the TL momenta K xy and 
K xz to zero and the TL momenta are reset to zero (Eqs. ()42)) - (j43|) ). Thus the remaining 
components contributing to weak hyperbolicity are controlled as a consequence of the mo- 
mentum constraint. This is of interest because the method of required the addition of 
the momentum constraint to control the behavior of the TL components in the evolution 
equations. 

VIII. RELATION TO STANDARD GAUGE SETTING FOR LINEARIZED 
WAVES 

In ("elementary") textbook formulations, it is common to say that a gauge can be set 
to put h yy + h zz = and h x i = initially, and that this gauge is then maintained by the 
evolution equations. Those formulations are set in terms of 4— dimensional gauge setting, 
rather than a 3 + 1 split as we employ. We have seen that h yy + h zz = is a requirement 
of the Hamiltonian constraint, modulo integration constants. In our 3 + 1 context, we then 
see that the h xx stays zero only if the initial slice is chosen K = 0. This is of course a choice 
of the initial slice, so this is a 4— dimensional gauge choice. 

This observation is relevant to the superficially very strange feature that, even just concen- 
trating on the h xx equation, strong hyperbolicity requires densitized lapse. Mathematically 
this is perfectly well defined and consistent. To a physicist, however, it is extremely strange 
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that we must warp our concept of t = constant to achieve controllable solutions, and that 
this process actually accomplishes mathematically well behaved solutions. However, in the 
context of the observation that K need not be initially set to zero, the situation becomes 
more comprehensible. The choice K ^ is not consistent with the elementary textbook 
gauge setting just described. Instead by taking K ^ one has taken a peculiarly curved 
initial spatial slice; the specific choice of densitized lapse is required to correct it. While it 
may still seem strange that a wavelike (propagating) solution is the best way to accomplish 
it, this is what happens. 

Further, another part of the assertion in elementary treatments is that the TL components 
h xy and h xz be set to zero as a gauge choice, and will remain zero. The approach from 
hyperbolicity achieves this by introducing another equation with contribution on the right 
side proportional to the momentum constraint, which is ideally zero. The resulting system 
maintains evolution close to satisfaction of the gauge conditions. The constrained evolution 
approach guarantees exact solution of this gauge condition (h xy and h xz remain constant, 
and can be set to zero by a simple rotation around the x— axis). Because small errors can 
arise in any computational system, a combination of the two approaches seems the best 
approach to achieving long term computational relativity code stability. 
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Appendix A: Generality of Study of Hyperbolicity 

The analysis of the hyperbolicity in section IV Bl is written in terms of linearized gravity 
as defined in section HTT1 However, based on arguments in Q|, the results are completely 
general for ADM, whet her p osed in first or second order form. We now summarize arguments 
for several points from 



her p c 
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(3 — 0: For reasonable shift vectors (shift velocities less than c), "the role of the shift 
vector is to displace the value of the eigenvalue... , and the lapse re-scales it. But a change 
of lapse... and shift can not change a real eigenvalue into an imaginary one. It can not affect 
the hyperbolicity of the system." See [l^ . section III A. 

Linearization: The hyperbolicity classification is based on the principal part of the sym- 
bol of the differential operator, which involves only the highest derivatives. Hence the terms 
shown in Eqs.(J7J), (JSJ), or equivalently Eq. ©, above define the hyperbolicity. This is true es- 
sentially because the Einstein equations are quasilinear, i.e. linear in the highest derivatives, 
and this is the reason that one can consider Fourier transforms of the linearized equations; 
consideration of hyperbolicity essentially looks at only the highest frequency behavior. 

We close a defect in the statement just made by following an argument given by Q], 
"since the norms | \$ and | \h are equivalent and smoothly related, ... the properties of the 
eigenvalues and eigenvectors of the principal symbol are the same with either norm" . Thus 
the undifferentiated metric (or inverse metric) coefficients in the Einstein equations can be 
equivalently replaced by Sij or its inverse, leading to Eqs. ((Zj), (jEJ) or the equivalent second 
order form Eq. (JHJ). 

Restriction to 1 — d Spatial Dependence: With the results just above, the background is 
isotropic, and nothing is lost in picking a specific direction for the wavevector hi associated 
with the local gradient. 

Note that we do not claim generality for application of constraint solution beyond the 
linearized wave spacetime. 



Appendix B: Other Boundary Conditions in Constraint Solution 

One of the advantages of writing explicit (if simple) solution of the linearized constraint 
equations is that one can explicitly investigate the effect of boundary conditions. For in- 
stance, in full 3 — d constraint solution we impose the inner conditions 0=1 and w l = 0, 
set at the excision surface of the black hole data. We take the outer conditions as in 
they are versions of the Robin ll| condition. Thus for the black hole case the boundary 
condition for <p is taken as d r (r(<j)— 1)) = 0, and the outer boundary condition for the radial 
component of w l , i.e. w r , is taken as d r (rw r ) = 0. Define also the transverse components 
of w l , X^ T = w l (5j — fjf- 7 ) (where the T means transverse). The outer boundary conditions 
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on X lT are d r (r 2 X lT ) = 0. We can model these boundary conditions in our 1 — d system by 
assuming a finite x— domain, say [0,x m ]. (As in the 3 — d situation, we will consider the 
outer boundary x m becoming large, x m — > oo.) We thus write the boundary conditions in 
terms of x— derivatives rather than r— derivatives. Since we want to impose conditions on 
w % rather than its derivatives, the solution for w l x , Eqs. - can be integrated once 
again by writing the explicit (0 to x) integration of the source terms, and introducing new 
integration constants Di. For instance 



w y 



f X A yx dx + C y x + D y (48) 
Jo 



It is straightforward to verify that one can consistently set each of the variables 0i, w l 
to zero at x = 0, and satisfy the Robin conditions at x = x m as we now show. (Notice that 
in our linearized problem, <ft — 1 = <pi.) 

We want w y = at x = and d x (x 2 w y ) = at x = x m . Thus D y = and 



2 r x m _ 

3C y = — / A yx dx + A yx {x m ). (49) 

Xyyi JO 

A completely analogous treatment applies to C z . 

For w x we take conditions w x = at x = (which gives D x = 0), and d x (xw x ) = at 
x = x m , which gives 

2C X = (A xx - \K){x m ) + — [ Xm (A xx - \K)dx. (50) 

The solution for <pi directly involves two integration constants, A, B, cf. Eq. ()30|) . We 
take inner boundary condition 0i = at x = 0, which sets B = \{h yy + h yy )\Q. We then set 
dx{.x4>i) — at x = x m , which gives 



16A = d x {h yy + h zz )(x m ) + —ihyy + ~h Z z)\o m - (51) 

•Em 

The ideal computational boundary condition has x m — > oo. With reasonable assumptions 
about the boundedness of the data at large x m , Eqs.fl4*9*]) - (joTIj) lead to Cj — > 
In Eq. (fSTjl . similarly A — > £lS ^ oo. 
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Thus the w l have precisely the form assumed in the main text with C\ — 0. For the 
situation is only slightly mor complicated. We have 

<Pi = -\{hyy + h zz )\l. (52) 
^From Eq. (JHHj) this gives the result for the metric components after the constraint solve: 

kyy + k ZZ = (hyy + H zz )q. (53) 

(54) 

Such a constant value can be removed by a gauge (coordinate) transformation transverse 
to the x propagation direction. Hence we may take (h yy + h zz ) = in setting the data. 

The result of the analysis in this Appendix is thus that the analysis in the main text 
accurately models the procedures used in our 3 — d simulations. This has the important 
corollary that the boundary conditions in the 3 — d case are completely internally consistent; 
we have not set too restrictive conditions (which might have only isolated solutions, i.e. 
might define an eigenvalue problem). 

It is not necessary to consider infinite domains x m — > oo. To treat cosmologies for instance, 
it is more appropriate to consider a specific finite value of x m . In this case the integration 
constants do not vanish. As an example, using Eq. (J25j) we find that the solution for A xx is 
A xx = %K + C X . Assuming that K changes little from timestep to timestep, the constrained 
evolution reproduces this same value of C x . Clearly C x can be consistently set to zero by 
setting A xx = \K exactly, regardless of the value of x m 

As above, we can interpret the integration constants from the Hamiltonian constraint as 
a kind of gauge choice in the initial data, which is re-imposed in the constrained evolution. 
For instance, the constants A, B, which arise from the Hamiltonian constraint, were set to 
zero in the body of the paper, implying that the transverse trace h yy + h zz is set to zero. The 
terms 8 Ax + 8B in Eq. (}35|) give constants to h yy + h zz , and the constant C x appearing in 
Eq. gives a time derivative to this transverse trace. Note that while h yy + h zz satisfies 
a wave equation, Eq. (fT3*|) . none of the integration constant terms has a second spatial 
derivative. Hence they will evolve trivially: 

h yy + h zz = 8Ax + 8B + 2C x t. (55) 
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With this form, it can be seen that repeated application of the constraints, as in constrained 
evolution, consistently leads to the same value of the constant A. The form Eq. (J55J) identifies 
the terms associated with the integration constants as gauge terms (the gauge vector £ x = 
2A(y 2 + z 2 ), ^ y = —AAxy, £ z = —4Axz removes the time independent terms). The same 
method can be used to remove the time dependent term at any one time, but this secular 
C x t term can eventually destroy the linearization assumption (unless C x is forced to zero 
as suggeted above). The secular term can be eliminated by an outgoing wave boundary 
condition, which forces C x = 0. 

Similar considerations apply to the terms C y and C z in Eqs. (|4"2|) . (|4Hj) . 

If the constraints are repeatedly applied, as in our form of constrained evolution, the 
value of C y is consistent. That is, from Eq. (J42j) . A xy = K xy has the constant value C y . The 
evolution equation Eq. (JTTJl shows that K xy remains constant (ideally) during an integration 
step. If K xy has the value C y as integrated in Eq. (|4"2j) . then Eq. (|4"5j) consisitently returns 
the same value of C y . Clearly the behavior of K xz is similar. Again, such terms can cause 
secular growth, but such growth can be controlled by, for instance, putting an outgoing wave 
boundary condition on h xy and h xz . 
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